# -*coding:utf-8 -*

import numpy as np
import matplotlib.pyplot as mp
from scipy.integrate import odeint

#------------------------#
#       projet 6         #
#                        #
#       partie 1         #
#------------------------#


# Implementation des differentes methodes de resolutions

# methode d'Euler
def step_euler(y,t,h,f) : 
    return  y + h * f(y,t)

#methode du point milieu
def step_midpoint(y,t,h,f):
    return y + h * f(y + (h * f(y,t) / 2), t+(h/2))

#methode de Heun
def step_heun(y,t,h,f) :
    return y + 0.5 * h * (f(y,t) + f(y+h*f(y,t) , t+h))

# methode de Runge-Kutta d'ordre 4    
def step_RK4(y,t,h,f) :
    k1 = f(y,t)
    k2 = f(y +(h * k1 / 2), t+(h/2))
    k3 = f( y + (h * k2 /2), t + (h/2))
    k4 = f( y+ h*k3, t+h)
    return y + h*(k1+ 2*k2 + 2*k3 + k4) /6.

#Premier algorithme de resolution
def meth_n_step(y0, t0, N,h,f, meth) :
    y = y0
    t = t0
    res = []
    for i in range (N) :
        res.append(y)
        y = meth(y,t,h,f)
        t += h
    return res


#Deuxieme algorithme de resolution
def meth_epsilon(y0, t0, tf, eps, f, meth):
    h = 1.0
    y2 = meth_n_step(y0, t0, np.rint(tf/h), h, f, meth)
    e = 10000000.0
    while (e > eps):
        y1 = y2
        h = h / 2
        y2 = meth_n_step(y0, t0, np.rint(tf/h), h, f, meth)
        max1 = len(y1)
        max2 = len(y2)
        e =y2[max2-1] - y1[max1-1]
    return y2


#Question 3
def champ_tangente(f,x1,x2,h,y0,meth):
    "fonction qui trace le champ des tangentes de l'équation différentielle passée en paramètre d'entrée"
    x = np.arange(x1,x2,h)
    y = np.arange(x1,x2,h)
    n = x.size
    
    for i in range (1,n) :
        for j in range(1,n) :
            X = np.array([[x[i]]])
            Y = np.array([[y[j]]])
            norme = np.sqrt(1.0 + (f(y[j],x[i]))**2)
            V = np.array([ [f(y[j],x[i])] /norme ])
            U = np.array([ [1.0] / norme ])
            mp.quiver(X,Y,U, V, color="k", width = .002, scale = 50)

    abs = np.arange(x1, x2, .01)
    ord = np.arange(x1, x2, .01)
    for i in range(0, len(ord)):
        abs[i] = 0.0
    mp.plot(abs, ord,color= "k", linewidth=1.0)

    t = np.arange(x1, x1 + len(y) * h, h)
   
    sol1 = meth_n_step(-0.5, x1, np.rint((x2-x1)/h), h, f, meth)
    curve1 =mp.plot (t, sol1, linewidth=1.0)

    sol2 = meth_n_step(0., x1, np.rint((x2-x1)/h), h, f, meth)
    curve2=mp.plot (t, sol2, linewidth=1.0)

    sol3 = meth_n_step(0.5, x1, np.rint((x2-x1)/h), h, f, meth)
    curve3=mp.plot (t, sol3, linewidth=1.0)
        
    mp.ylim([x1,x2])
    mp.xlim([x1,x2])
    mp.xlabel('t')
    mp.ylabel('y')
    mp.title("Graphe du champ des tangentes de l'equation differentielle")
    mp.legend((curve1,curve2,curve3),('y0 = -0.5','y0 = 0','y0 = 0.5'))
    mp.show()

# Question 4
            

#fonction qui trace les courbes obtenues et les comprae avec la solution theorique
def graphe_dim1(f,y0,t0,h):
    mp.title("Tracage des courbes de l'equation differentielle")
    mp.xlabel("t")
    mp.ylabel("y")
  
    de = meth_n_step(y0,t0,100,h,f,step_euler)
    dm = meth_n_step(y0,t0,100,h,f,step_midpoint)
    dh = meth_n_step(y0,t0,100,h,f,step_heun)
    d4 = meth_n_step(y0,t0,100,h,f,step_RK4)
    t = np.linspace(t0,10.,len(de))
    dt = odeint(f,1.,t)    
    
    e_curve = mp.plot(t,de)
    m_curve = mp.plot(t,dm)
    h_curve = mp.plot(t,dh)
    r_curve = mp.plot(t,d4)
    t_curve = mp.plot(t,dt)
    
    mp.legend((e_curve,m_curve,h_curve,r_curve, t_curve),('Euler','point milieu','Heun', 'Runge Kutta 4','solution theorique'))
    mp.show()


#============= fonction qui trace la solution en dimension 2 =======#

def graphe_dim2 (f,y0,t0,h):
    mp.title("Tracage des courbes de l'equation differentielle en dimension 2")
    mp.xlabel("t")
    mp.ylabel("y")
  
    de = meth_n_step(y0,t0,100,h,f,step_euler)
    dm = meth_n_step(y0,t0,100,h,f,step_midpoint)
    dh = meth_n_step(y0,t0,100,h,f,step_heun)
    d4 = meth_n_step(y0,t0,100,h,f,step_RK4)

    t = np.linspace(t0,10.,len(de))
    
    #e_curve = mp.plot(t,de)
    #m_curve = mp.plot(t,dm)
    h_curve = mp.plot(t,dh)
    r_curve = mp.plot(t,d4)
    
    mp.legend((h_curve,r_curve),('Heun (1)', 'Heun (2)','Runge Kutta 4 (1)','Runge Kutta (2)'))
    mp.show()


#=============fonction qui trace y2 en fonction de y1============#

def graphe2(f,h):

    y1_euler = np.arange(0.0, 10.0, h)
    y2_euler = np.arange(0.0, 10.0, h)
    y1_milieu = np.arange(0.0, 10.0, h)
    y2_milieu = np.arange(0.0, 10.0, h)
    y1_heun = np.arange(0.0, 10.0, h)
    y2_heun = np.arange(0.0, 10.0, h)
    y1_runge = np.arange(0.0, 10.0, h)
    y2_runge = np.arange(0.0, 10.0, h)

    

    y1_euler[0] = np.array([1.0]) 
    y2_euler[0] = np.array([0.0])
    y1_milieu[0] = np.array([1.0]) 
    y2_milieu[0] = np.array([0.0])
    y1_heun[0] = np.array([1.0]) 
    y2_heun[0] = np.array([0.0])
    y1_runge[0] = np.array([1.0]) 
    y2_runge[0] = np.array([0.0])
    

    for i in range (1, len(y1_euler)):
        y_euler = np.array([y1_euler[i-1], y2_euler[i-1]])
        vect_euler = step_euler(y_euler, 0.0, h, f)
        [y1_euler[i], y2_euler[i]] = vect_euler


        y_milieu = np.array([y1_milieu[i-1], y2_milieu[i-1]])
        vect_milieu = step_midpoint(y_milieu, 0.0, h, f)
        [y1_milieu[i], y2_milieu[i]] = vect_milieu


        y_heun = np.array([y1_heun[i-1], y2_heun[i-1]])
        vect_heun = step_heun(y_heun, 0.0, h, f)
        [y1_heun[i], y2_heun[i]] = vect_heun


        y_runge = np.array([y1_runge[i-1], y2_runge[i-1]])
        vect_runge = step_RK4(y_runge, 0.0, h, f)
        [y1_runge[i], y2_runge[i]] = vect_runge


    mp.plot (y1_euler, y2_euler, 'r', label = 'Euler')
    mp.plot (y1_milieu, y2_milieu, 'b', label = 'Point milieu')
    mp.plot (y1_heun, y2_heun, 'g', label = 'Heun')
    mp.plot (y1_runge, y2_runge, 'k', label = 'Runge-Kutta')



    y1 = np.arange(0.0, 10.0, .01)
    y2 = np.arange(0.0, 10.0, .01)
    t = np.arange(0.0, 10.0, .01)
    for i in range (0, len(t)):
        y1[i] = np.cos(t[i])
        y2[i] = np.sin(t[i])
    mp.plot (y1, y2, label = 'Solution exacte')

    mp.xlabel('y1')
    mp.ylabel('y2')
    mp.legend()
    mp.title("y2 en fonction de y1")
    mp.show()


